## AR(3) model with shift in mean

rm(list=ls())
set.seed(616)
graphics.off()
library(sandwich)
library(dplyr)
source("R/data_fns.R")
source("R/sdr_fns.R")

########################################
## CHOOSE HERE

spec <- "y1"
data <- loadShortRate()

## spec <- "y10"
## data <- loadLongRate()

########################################

byear <- 1990           # break year
N <- 400                # maximum maturity
mats <- 1:N             # yield maturities in years (for plotting)

## subtract out different unconditional mean pre and post break
d1 <- data %>% filter(year <= byear)
d2 <- data %>% filter(year > byear)
(mean1 <- d1 %>% pull(rr) %>% mean)
(mean2 <- d2 %>% pull(rr) %>% mean)
data <- mutate(data,
               rtilde = rr - ifelse(year <= byear, mean1, mean2))
mm <- lm(rr ~ 1, d1)
cat("Mean and SE before break:", etabar1 <- mm$coef, etase1 <- sqrt(NeweyWest(mm))[1], "\n")
mm <- lm(rr ~ 1, d2)
cat("Mean and SE after break:", etabar2 <- mm$coef, etase2 <- sqrt(NeweyWest(mm))[1], "\n")
stopifnot(all.equal(mean1, etabar1, check.attributes=FALSE))
stopifnot(all.equal(mean2, etabar2, check.attributes=FALSE))
mm <- lm(rr ~ I(year <= byear), data)
stopifnot(all.equal(mean2, mm$coef[1], check.attributes=FALSE))
cat("t-stat on mean change:", mm$coef[2]/sqrt(diag(NeweyWest(mm)))[2], "\n")

## AR model on residuals
mod <- ar.ols(data$rtilde, aic=FALSE, order.max=3, demean=FALSE)
rhobar <- mod$ar[,,]
rhose <- mod$asy.se.coef$ar
cat("# Sum of AR coeffs:", sum(rhobar), "\n")
sig <- sqrt(mod$var.pred)
tbl <- cbind(mod$ar[,,], mod$asy.se.coef$ar)
tbl <- cbind(tbl, tbl[,1]/tbl[,2])
colnames(tbl) <- c("Estimate", "SE", "t-stat")
rownames(tbl) <- paste("lag", 1:mod$order)
print(round(tbl, 5))

## simulate SDR to get term structure and save results
P1 <- sim_ar_sdr(etabar1, etase1, rhobar, rhose, sig)
P2 <- sim_ar_sdr(etabar2, etase2, rhobar, rhose, sig)
y1 <- -100/mats*log(P1)
y2 <- -100/mats*log(P2)
rstar_old <- mean1
rstar_new <- mean2
save(file=paste0("results/sdr_meanshift_",spec,".RData"), rstar_old, rstar_new, P1, y1, P2, y2) #

dev.new()
par(mar=c(3,3,2,.5), mgp=c(2,.6,0))
plot(mats, y1, type="l", ylim=range(0, y1, y2, 5), lwd=2, ylab="Discount rate", xlab="Horizon (years)", col="red", main="SDRs from AR(3) with breaking mean")
lines(mats, y2, lwd=2, col="steelblue")
abline(h=0, lty=2)
abline(h=mean1, lty=2, col="red", lwd=2)
abline(h=mean2, lty=2, col="steelblue", lwd=2)
legend("topright", c("pre-1990 mean", "post-1990 mean"), col=c("red", "steelblue"), lwd=2, bg="white")

